%% Preprocess and choose initial centers of holes

%x1=read image you want to process
N=4;
h=zeros(2*N+1);
for i=1:2*N+1
    for j=1:2*N+1
       h(i,j)=mvnpdf([i-N-1 j-N-1],[0 0],[1 1]); 
   
    end
end
h=h/sum(sum(h));  
y1=conv2(x1,h,'same'); %low pass filtering

% determine initial centers 
y3=imbinarize(y1,'adaptive'); %or (y1>threshold);
se=strel('disk',2);
y4=imopen(y3,se);
figure;imshow(y4,[]);

figure;imagesc(x1);hold on;
[B,L] = bwboundaries(y4,"noholes");
stats = regionprops(L,"Circularity","Centroid","MajorAxisLength","MinorAxisLength");
cendoto=zeros(length(B),2);
radii=mean([stats.MajorAxisLength stats.MinorAxisLength],2)/2;
 
for k = 1:length(B)
  boundary = B{k};
  centroid = stats(k).Centroid;
  cendoto(k,:)=centroid;
  plot(centroid(1),centroid(2),"r.",'MarkerSize',15);
end
defects=[];

%% Add and delete some centers manually
figure;imagesc(x1);hold on;
scatter(cendoto(:,1),cendoto(:,2),10,'ro','filled');
hold on;
while strcmp(questdlg('Do you want to delete some points?'),'Yes')
    [x,y]=ginput;
    [~,pn]=min((cendoto(:,1)-x).^2+(cendoto(:,2)-y).^2);
    plot(cendoto(pn,1),cendoto(pn,2),"c.",'MarkerSize',15);
    cendoto(pn,:)=[];
end
while strcmp(questdlg('Do you want to add some points?'),'Yes')
    [x,y]=ginput;
     plot(x,y,"r.",'MarkerSize',15);
     cendoto=[cendoto;x,y];
end

while strcmp(questdlg('Do you want to add some black defects?'),'Yes')
    [x,y]=ginput;
     plot(x,y,"w.",'MarkerSize',15);
     defects=[defects;x,y];
end
inicendot=[cendoto;defects];

%% iteration for determining the final centers of holes
A1=(1:size(x1,2))';
B1=ones(size(x1,1),1);
C=zeros(size(x1,1)*size(x1,2),2);
C(:,1)=kron(A1,B1);
C(:,2)=repmat((1:size(x1,1))',size(x1,2),1);
clasf=zeros(size(x1,1)*size(x1,2),1);
Change=0;
Changeline=[];
cthreshold=0.1;
cendot=inicendot;
for k=1:80 % iterate for 80 steps
    for i=1:size(x1,1)*size(x1,2)
        dis=(C(i,1)-cendot(:,1)).^2+(C(i,2)-cendot(:,2)).^2;
        [sdis,I]=sort(dis);
        clasf(i,1)=I(1);
    end
    
    ncendot=zeros(size(cendot,1),2);
    for i=1:size(cendot,1)
        aroundot=C(clasf==i,:);
        if i<size(cendoto,1) % for bright holes, find the centers of bright regions
            val=(x1(sub2ind(size(x1),aroundot(:,2),aroundot(:,1))));
        else % for dark defects,find the centers of dark regions
            val=max(max(x1))-(x1(sub2ind(size(x1),aroundot(:,2),aroundot(:,1))));
        end
        w=val;
        ncendot(i,1)=sum(aroundot(:,1).*w)/sum(w);
        ncendot(i,2)=sum(aroundot(:,2).*w)/sum(w);
    end  
    Change=sum(sqrt((ncendot(:,1)-cendot(:,1)).^2+(ncendot(:,2)-cendot(:,2)).^2));
    Changeline=[Changeline,Change];  %record changes between every two steps
    cendot=ncendot;
    if Change<cthreshold %threshold for ending the iterations
        break;
    end
end
figure;plot(Changeline);
figure;imshow(x1);hold on;
for i=1:size(cendot(:,1))
  scatter(cendot(i,1),cendot(i,2),45,'ro','filled');
end
voronoi(cendot(:,1),cendot(:,2)); % draw the Voronoi cells

%% count nearest points (to draw the hole lattice in Fig.S2E)
ncendot=cendot;
figure;imshow(x1,[]);hold on;[vx,vy]=voronoi(ncendot(:,1),ncendot(:,2));

crossp1=[vx(1,:)',vy(1,:)'];
crossp2=[vx(2,:)',vy(2,:)'];
[crossp,~,~] = union(crossp1,crossp2,'rows');
crossp=crossp.*[(crossp(:,1)>0).*(crossp(:,1)<size(x1,2)),(crossp(:,1)>0).*(crossp(:,1)<size(x1,2))];
crossp=crossp.*[(crossp(:,2)>0).*(crossp(:,2)<size(x1,1)),(crossp(:,2)>0).*(crossp(:,2)<size(x1,1))];
crossp(any(crossp,2)==0,:)=[];
pointcount=zeros(size(cendot,1),1);

nearpoint=cell(size(cendot,1),1);
for i=1:size(crossp,1)
    dis2=sqrt((crossp(i,1)-cendot(:,1)).^2+(crossp(i,2)-cendot(:,2)).^2);
    kk=sum(round(dis2,2)==round(min(dis2),2));
    pointcount=pointcount+(round(dis2,2)==round(min(dis2),2));
    pn=(find((round(dis2,2)==round(min(dis2),2))));
    for j=1:size(pn)
       nearpoint(pn(j),1)={unique([cell2mat(nearpoint(pn(j),1)),pn'])};
    end

end

%% calculate the local density of holes (for figure 4)
areasize=zeros(size(cendot,1),1);
colorm=zeros(size(x1,1),size(x1,2));
for i=1:size(x1,1)
    for j=1:size(x1,2)
        [~,pn]=min((cendot(:,1)-j).^2+(cendot(:,2)-i).^2);        
        areasize(pn)=areasize(pn)+1;
    end
end
for i=1:size(x1,1)
    for j=1:size(x1,2)
        [~,pn]=min((cendot(:,1)-j).^2+(cendot(:,2)-i).^2);        
        colorm(i,j)=areasize(pn);
    end
end
%asize= real size for each pixel in the unit of nm
colorm=1e14./(colorm*asize);
